Fit stomach content models

Load libraries

# Load libraries, install if needed
pkgs <- c("tidyverse", "tidylog", "devtools", "sdmTMB", "sdmTMBextra", "terra", "mapplots",
          "viridis", "visreg", "modelr", "future", "kableExtra", "ggh4x", "patchwork") 

if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library, character.only = T))

# Import some plotting functions
# Source code for map plots
# You need: # devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/pred-prey-overlap/main/R/functions/map-plot.R")
options(ggplot2.continuous.colour = "viridis")

# devtools::install_github("seananderson/ggsidekick") # not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Set path
home <- here::here()

# For crossvalidation: paralell processing
plan(multisession)

set.seed(99) 

# To load entire cache in interactive r session, do: 
# qwraps2::lazyload_cache_dir(path = paste0(home, "/R/main-analysis/03-fit-diet-models_cache/html"))

Read cleaned data

df <- read_csv(paste0(home, "/data/clean/stomachs.csv")) |>
  #filter(pred_length > 15 & pred_length < 50) |> 
  mutate(depth_sc = (depth - mean(depth))/sd(depth),
         year_f = as.factor(year),
         month_f = as.factor(month),
         ices_rect = as.factor(ices_rect),
         pred_length_sc = (pred_length - mean(pred_length)) / sd(pred_length)) # read_csv(paste0(home, "/data/clean/stomachs.csv")) |> summarise(mean = mean(pred_length))

# Read the prediction grid...
pred_grid <- bind_rows(read_csv(paste0(home, "/data/clean/pred_grid_(1_2).csv")),
                       read_csv(paste0(home, "/data/clean/pred_grid_(2_2).csv")))

# Summarize samples size per location for plotting... 
dd_plot <- df |>
  group_by(year, X, Y) |> 
  mutate(sample_size = n(),
         pos_id = paste(year, X, Y)) |> 
  ungroup() |> 
  distinct(pos_id, .keep_all = TRUE)

# Scale with respect to data!
df |> group_by(month) |> summarise(n = n()) |> arrange(desc(n))

# A tibble: 8 × 2
  month     n
  <dbl> <int>
1     3  7040
2    11  2207
3    12  1980
4     4   553
5     2   549
6    10    57
7     6    36
8     5     3

pred_grid <- pred_grid |> 
  drop_na(depth) |> 
  filter(depth < quantile(df$depth, probs = 0.99)) |> 
  filter(quarter == 4) |> # Not needed in theory for saduria...
  mutate(depth_sc = (depth - mean(df$depth)) / sd(df$depth),
         year_f = as.factor(year),
         month_f = as.factor(11),
         year = as.integer(year),
         pred_length_sc = 0,
         ices_rect = as.factor(ices_rect)) |> 
  droplevels()

ggplot() + 
  geom_histogram(data = pred_grid, aes(depth, fill = "pred_grid")) + 
  geom_histogram(data = df, aes(depth, fill = "data"))

ggplot(pred_grid |> filter(year == 1999), aes(X, Y, color = depth)) + 
  geom_point()

Plot data!

df |> 
  pivot_longer(c(FR_spr, FR_her, FR_sad)) |> 
  ggplot(aes(year, value)) +
  geom_jitter(height = 0, alpha = 0.5) + 
  coord_cartesian(ylim = c(0, 0.1)) +
  facet_wrap(~name, ncol = 1) +
  geom_smooth()

Set up a mesh

mesh <- make_mesh(df, c("X", "Y"), cutoff = 6)

ggplot() +
  inlabru::gg(mesh$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = df, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.3, vjust = 3) + 
  labs(x = "Easting (km)", y = "Northing (km)")

Spatial cross validation and AIC to compare spatial vs non spatial models

5-fold spatial cross validation & AIC to compare two models: spatial fields or ices rectangle as random effects

# Set up spatial clusters
clust <- kmeans(df[, c("X", "Y")], 5)$cluster

# Plot
df$clust_id <- clust

plot_map +
  geom_point(data = df, aes(X*1000, Y*1000, color = as.factor(clust_id))) +
  scale_color_viridis(discrete = TRUE) + 
  geom_sf(size = 0.1)

Warning: Removed 39 rows containing missing values (`geom_point()`).

Warning: colourbar guide needs continuous scales.

# Sprat
# ICES as a random effects
spr_space_1_cv <- sdmTMB_cv(
  FR_spr ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
  data = df,
  mesh = mesh,
  spatial = "off",
  family = tweedie(link = "log"),
  fold_ids = clust,
  k_folds = length(unique(clust))
)

Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.

# Spatial random field
spr_space_2_cv <- sdmTMB_cv(
  FR_spr ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
  data = df,
  mesh = mesh,
  spatial = "on",
  family = tweedie(link = "log"),
  fold_ids = clust,
  k_folds = length(unique(clust))
)

Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.

# Herring
# ICES as a random effects
her_space_1_cv <- sdmTMB_cv(
  FR_her ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
  data = df,
  mesh = mesh,
  spatial = "off",
  family = tweedie(link = "log"),
  fold_ids = clust,
  k_folds = length(unique(clust))
)

Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.

# Spatial random field
her_space_2_cv <- sdmTMB_cv(
  FR_her ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
  data = df,
  mesh = mesh,
  spatial = "on",
  family = tweedie(link = "log"),
  fold_ids = clust,
  k_folds = length(unique(clust))
)

Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.

# Saduria
# ICES as a random effects
sad_space_1_cv <- sdmTMB_cv(
  FR_sad ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect),
  data = df,
  mesh = mesh,
  spatial = "off",
  family = tweedie(link = "log"),
  fold_ids = clust,
  k_folds = length(unique(clust))
)

Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.

# Spatial random field
sad_space_2_cv <- sdmTMB_cv(
  FR_sad ~ 0 + year_f + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
  data = df,
  mesh = mesh,
  spatial = "on",
  family = tweedie(link = "log"),
  fold_ids = clust,
  k_folds = length(unique(clust))
)

Running fits with `future.apply()`.
Set a parallel `future::plan()` to use parallel processing.

AIC

# ICES as a random effects
fit_spr_m1 <- sdmTMB(
  FR_spr ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect), 
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  spatial = "off",
  spatiotemporal = "off",
  family = tweedie(link = "log"))

Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`

sanity(fit_spr_m1)

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100

fit_her_m1 <- sdmTMB(
  FR_her ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect), 
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  spatial = "off",
  spatiotemporal = "off",
  family = tweedie(link = "log"))

Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`

sanity(fit_her_m1)

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100

fit_sad_m1 <- sdmTMB(
  FR_sad ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f) + (1|ices_rect), 
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  spatial = "off",
  spatiotemporal = "off",
  family = tweedie(link = "log"))

Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`

sanity(fit_sad_m1)

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✖ `ln_smooth_sigma` standard error may be large
ℹ Try simplifying the model, adjusting the mesh, or adding priors

✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100

# Spatial random field
fit_spr_m2 <- sdmTMB(
  FR_spr ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f), 
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  spatial = "on",
  spatiotemporal = "off",
  family = tweedie(link = "log"))

Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`

sanity(fit_spr_m2)

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large

fit_her_m2 <- sdmTMB(
  FR_her ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f), 
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  spatial = "on",
  spatiotemporal = "off",
  family = tweedie(link = "log"))

Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`

sanity(fit_her_m2)

✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large

fit_sad_m2 <- sdmTMB(
  FR_sad ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f), 
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  spatial = "on",
  spatiotemporal = "off",
  family = tweedie(link = "log"))

Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`

Print AIC & cross validation table

# Higher is better!
# likelihood from spatial cv
sum_loglik_space <- data.frame(Prey = c("Sprat", "Herring", "Saduria"),
                               rec  =  c(spr_space_1_cv$sum_loglik / 10000000,
                                         her_space_1_cv$sum_loglik / 10000000,
                                         sad_space_1_cv$sum_loglik / 10000000),
                               space = c(spr_space_2_cv$sum_loglik / 10000000,
                                         her_space_2_cv$sum_loglik / 10000000,
                                         sad_space_2_cv$sum_loglik / 10000000)) |> 
  mutate(rec_temp = rec, space_temp = space,
         rec = ifelse(rec_temp > space_temp, paste0("**", rec, "**"), rec),
         space = ifelse(space_temp > rec_temp, paste0("**", space, "**"), space)) |> 
  rename("sum loglik<sub>ICES rect</sub>" = rec,
         "sum loglik<sub>spatial</sub>" = space) |> 
  dplyr::select(-rec_temp, -space_temp)

sum_loglik_space

     Prey sum loglik<sub>ICES rect</sub> sum loglik<sub>spatial</sub>
1   Sprat       **5.58014982903663e-05**         4.90488173952715e-05
2 Herring             -0.850543167016049       **-0.807858283667988**
3 Saduria              -13.6423530384996        **-9.92716500451823**

# AIC ices models
aic_ices <- AIC(fit_spr_m1, fit_her_m1, fit_sad_m1) |>
  rownames_to_column() |> 
  mutate(Prey = "Sprat",
         Prey = ifelse(rowname == "fit_her_m1", "Herring", Prey),
         Prey = ifelse(rowname == "fit_sad_m1", "Saduria", Prey)) |> 
  mutate(AIC_r = AIC, 
         AIC_r2 = AIC) |> 
  dplyr::select(-rowname, -AIC, -df)

aic_space <- AIC(fit_spr_m2, fit_her_m2, fit_sad_m2) |>
  rownames_to_column() |> 
  mutate(Prey = "Sprat",
         Prey = ifelse(rowname == "fit_her_m2", "Herring", Prey),
         Prey = ifelse(rowname == "fit_sad_m2", "Saduria", Prey)) |> 
  
  mutate(AIC_s = AIC,
         AIC_s2 = AIC) |> 
  dplyr::select(-rowname, -AIC, -df)

aic <- aic_ices |> 
  left_join(aic_space, by = "Prey") |> 
  mutate(AIC_r = ifelse(AIC_r2 < AIC_s2, paste0("**", AIC_r, "**"), AIC_r),
         AIC_s = ifelse(AIC_s2 < AIC_r2, paste0("**", AIC_s, "**"), AIC_s)) |> 
  rename("AIC<sub>spatial</sub>" = AIC_s,
         "AIC<sub>ICES rect</sub>" = AIC_r) |> 
  dplyr::select(-AIC_s2, -AIC_r2)

aic

     Prey AIC<sub>ICES rect</sub> AIC<sub>spatial</sub>
1   Sprat              -2925.7767 **-3284.89068021079**
2 Herring                558.8468  **426.356157207818**
3 Saduria              -1323.1981 **-1547.59615868862**

kableExtra::kbl(sum_loglik_space, format = "markdown", caption = "5-fold Spatial cross validation", digits = 2) |>
  kable_styling(font_size = 20)

5-fold Spatial cross validation

Prey	sum loglikICES rect	sum loglikspatial
Sprat	5.58014982903663e-05	4.90488173952715e-05
Herring	-0.850543167016049	-0.807858283667988
Saduria	-13.6423530384996	-9.92716500451823

kableExtra::kbl(aic, format = "markdown", caption = "AIC", digits = 2) |>
  kable_styling(font_size = 20)

AIC

Prey	AICICES rect	AICspatial
Sprat	-2925.78	-3284.89068021079
Herring	558.85	426.356157207818
Saduria	-1323.20	-1547.59615868862

# TODO: it seems perhaps ?kbl doesn't obey rounding when I use **. Make table in word instead for the report

Check residuals and from selected models

# Summary
summary(fit_spr_m2)

Spatial model fit by ML ['sdmTMB']
Formula: FR_spr ~ 0 + depth_sc + s(pred_length_sc, k = 3) + (1 | month_f)
Mesh: mesh (isotropic covariance)
Time column: year
Data: df
Family: tweedie(link = 'log')
 
                coef.est coef.se
depth_sc            0.37    0.07
spred_length_sc     0.52    0.07

Smooth terms:
                    Std. Dev.
sds(pred_length_sc)     35.52

Random intercepts:
        Std. Dev.
month_f      0.82

Time-varying parameters:
                 coef.est coef.se
(Intercept)-1993    -6.45    0.43
(Intercept)-1994    -6.22    0.46
(Intercept)-1995    -7.38    0.43
(Intercept)-1996    -6.08    0.39
(Intercept)-1997    -7.31    0.45
(Intercept)-1998    -5.85    0.49
(Intercept)-1999    -6.20    0.38
(Intercept)-2000    -6.27    0.39
(Intercept)-2001    -5.66    0.40
(Intercept)-2002    -7.70    0.46
(Intercept)-2003    -4.69    0.40
(Intercept)-2004    -7.35    0.39
(Intercept)-2005    -6.20    0.38
(Intercept)-2006    -6.14    0.38
(Intercept)-2007    -6.35    0.37
(Intercept)-2008    -6.25    0.37
(Intercept)-2009    -5.90    0.38
(Intercept)-2012    -4.94    0.38
(Intercept)-2013    -5.62    0.37
(Intercept)-2014    -5.16    0.42
(Intercept)-2018    -6.00    0.39
(Intercept)-2021    -5.88    0.37

Dispersion parameter: 0.30
Tweedie p: 1.44
Matérn range: 16.14
Spatial SD: 1.24
ML criterion at convergence: -1651.445

See ?tidy.sdmTMB to extract these values as a data frame.

summary(fit_her_m2)

Spatial model fit by ML ['sdmTMB']
Formula: FR_her ~ 0 + depth_sc + s(pred_length_sc, k = 3) + (1 | month_f)
Mesh: mesh (isotropic covariance)
Time column: year
Data: df
Family: tweedie(link = 'log')
 
                coef.est coef.se
depth_sc           -0.23    0.09
spred_length_sc    -0.79    0.07

Smooth terms:
                    Std. Dev.
sds(pred_length_sc)      34.7

Random intercepts:
        Std. Dev.
month_f      0.12

Time-varying parameters:
                 coef.est coef.se
(Intercept)-1993    -9.63    0.49
(Intercept)-1994    -8.07    0.45
(Intercept)-1995    -7.29    0.37
(Intercept)-1996    -8.16    0.37
(Intercept)-1997    -6.92    0.35
(Intercept)-1998    -6.93    0.46
(Intercept)-1999    -7.62    0.31
(Intercept)-2000    -7.49    0.33
(Intercept)-2001    -7.40    0.35
(Intercept)-2002    -8.26    0.47
(Intercept)-2003    -5.66    0.28
(Intercept)-2004    -7.27    0.28
(Intercept)-2005    -8.15    0.31
(Intercept)-2006    -7.49    0.27
(Intercept)-2007    -7.77    0.28
(Intercept)-2008    -7.48    0.26
(Intercept)-2009    -6.65    0.29
(Intercept)-2012    -6.35    0.32
(Intercept)-2013    -6.29    0.28
(Intercept)-2014    -7.68    0.49
(Intercept)-2018    -8.13    0.32
(Intercept)-2021    -8.17    0.25

Dispersion parameter: 0.61
Tweedie p: 1.48
Matérn range: 16.11
Spatial SD: 1.28
ML criterion at convergence: 204.178

See ?tidy.sdmTMB to extract these values as a data frame.

summary(fit_sad_m2)

Spatial model fit by ML ['sdmTMB']
Formula: FR_sad ~ 0 + depth_sc + s(pred_length_sc, k = 3) + (1 | month_f)
Mesh: mesh (isotropic covariance)
Time column: year
Data: df
Family: tweedie(link = 'log')
 
                coef.est coef.se
depth_sc           -0.01    0.15
spred_length_sc     0.19    0.05

Smooth terms:
                    Std. Dev.
sds(pred_length_sc)         0

Random intercepts:
        Std. Dev.
month_f      0.95

Time-varying parameters:
                 coef.est coef.se
(Intercept)-1993    -9.19    1.08
(Intercept)-1994    -9.22    1.08
(Intercept)-1995    -8.63    1.01
(Intercept)-1996    -8.61    1.00
(Intercept)-1997    -8.61    1.01
(Intercept)-1998    -9.71    1.11
(Intercept)-1999    -8.58    0.99
(Intercept)-2000    -8.16    0.98
(Intercept)-2001    -7.53    0.98
(Intercept)-2002    -9.07    1.00
(Intercept)-2003    -8.43    0.98
(Intercept)-2004    -8.44    0.98
(Intercept)-2005    -8.24    0.98
(Intercept)-2006    -8.01    0.97
(Intercept)-2007    -8.80    0.98
(Intercept)-2008    -8.67    0.98
(Intercept)-2009    -8.61    0.99
(Intercept)-2012    -7.96    0.99
(Intercept)-2013    -8.48    0.99
(Intercept)-2014    -8.22    0.98
(Intercept)-2018    -7.54    0.95
(Intercept)-2021    -7.14    0.94

Dispersion parameter: 0.73
Tweedie p: 1.54
Matérn range: 80.79
Spatial SD: 2.99
ML criterion at convergence: -782.798

See ?tidy.sdmTMB to extract these values as a data frame.

**Possible issues detected! Check output of sanity().**

# Residuals
spr_res <- mcmc_res <- residuals(fit_spr_m2, type = "mle-mcmc", 
                                 mcmc_samples = sdmTMBextra::predict_mle_mcmc(fit_spr_m2,
                                                                              mcmc_iter = 201,
                                                                              mcmc_warmup = 200))

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.002765 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 27.65 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 32.4323 seconds (Warm-up)
Chain 1:                0.064402 seconds (Sampling)
Chain 1:                32.4967 seconds (Total)
Chain 1: 

df$spr_res <- as.vector(spr_res)

her_res <- mcmc_res <- residuals(fit_her_m2, type = "mle-mcmc",
                                 mcmc_samples = sdmTMBextra::predict_mle_mcmc(fit_her_m2,
                                                                              mcmc_iter = 201,
                                                                              mcmc_warmup = 200))

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.002747 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 27.47 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 19.0425 seconds (Warm-up)
Chain 1:                0.030647 seconds (Sampling)
Chain 1:                19.0732 seconds (Total)
Chain 1: 

df$her_res <- as.vector(her_res)

sad_res <- mcmc_res <- residuals(fit_sad_m2, type = "mle-mcmc",
                                 mcmc_samples = sdmTMBextra::predict_mle_mcmc(fit_sad_m2,
                                                                              mcmc_iter = 201,
                                                                              mcmc_warmup = 200))

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.00279 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 27.9 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 243.25 seconds (Warm-up)
Chain 1:                0.251266 seconds (Sampling)
Chain 1:                243.501 seconds (Total)
Chain 1: 

df$sad_res <- as.vector(sad_res)

# Plot all together
df |> 
  rename(Sprat = spr_res,
         Herring = her_res,
         Saduria = sad_res) |> 
  pivot_longer(c(Sprat, Herring, Saduria)) |> 
  ggplot(aes(sample = value)) +
  facet_wrap(~name) +
  stat_qq(size = 0.75, shape = 21, fill = NA) +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
  theme(aspect.ratio = 1)

rename: renamed 3 variables (Sprat, Herring, Saduria)

pivot_longer: reorganized (Sprat, Herring, Saduria) into (name, value) [was 12425x25, now 37275x24]

ggsave(paste0(home, "/figures/supp/diet_qq.pdf"), width = 17, height = 11, units = "cm")

Plot conditional effects, random effects and spatial predictions

# Depth
vis_spr_dep <- visreg(fit_spr_m2, xvar = "depth_sc", plot = FALSE)
vis_her_dep <- visreg(fit_her_m2, xvar = "depth_sc", plot = FALSE)
vis_sad_dep <- visreg(fit_sad_m2, xvar = "depth_sc", plot = FALSE)

vis_dep <- bind_rows(vis_spr_dep$fit |> mutate(species = "Sprat"),
                     vis_her_dep$fit |> mutate(species = "Herring"), 
                     vis_sad_dep$fit |> mutate(species = "Saduria")) |> 
  mutate(var = "Depth (scaled)") |> 
  rename(x = depth_sc)

d_dep <- bind_rows(vis_spr_dep$res |> mutate(species = "Sprat"),
                   vis_her_dep$res |> mutate(species = "Herring"),
                   vis_sad_dep$res |> mutate(species = "Saduria")) |> 
  mutate(var = "Depth (scaled)") |> 
  rename(x = depth_sc)

# Month
vis_spr_mon <- visreg(fit_spr_m2, xvar = "month_f", plot = FALSE)
vis_her_mon <- visreg(fit_her_m2, xvar = "month_f", plot = FALSE)
vis_sad_mon <- visreg(fit_sad_m2, xvar = "month_f", plot = FALSE)

vis_mon <- bind_rows(vis_spr_mon$fit |> mutate(species = "Sprat"),
                     vis_her_mon$fit |> mutate(species = "Herring"), 
                     vis_sad_mon$fit |> mutate(species = "Saduria")) |> 
  mutate(var = "Month") |> 
  rename(x = month_f) |> 
  mutate(x = as.numeric(as.character(x)))

d_mon <- bind_rows(vis_spr_mon$res |> mutate(species = "Sprat"),
                   vis_her_mon$res |> mutate(species = "Herring"),
                   vis_sad_mon$res |> mutate(species = "Saduria")) |> 
  mutate(var = "Month") |> 
  rename(x = month_f) |> 
  mutate(x = as.numeric(as.character(x)))

# Predator length
vis_spr_len <- visreg(fit_spr_m2, xvar = "pred_length_sc", plot = FALSE)
vis_her_len <- visreg(fit_her_m2, xvar = "pred_length_sc", plot = FALSE)
vis_sad_len <- visreg(fit_sad_m2, xvar = "pred_length_sc", plot = FALSE)

vis_len <- bind_rows(vis_spr_len$fit |> mutate(species = "Sprat"),
                     vis_her_len$fit |> mutate(species = "Herring"),
                     vis_sad_len$fit |> mutate(species = "Saduria")) |> 
  mutate(var = "Predator length") |> 
  rename(x = pred_length_sc)

d_len <- bind_rows(vis_spr_len$res |> mutate(species = "Sprat"),
                   vis_her_len$res |> mutate(species = "Herring"),
                   vis_sad_len$res |> mutate(species = "Saduria")) |> 
  mutate(var = "Predator length") |> 
  rename(x = pred_length_sc)

vis <- bind_rows(vis_dep, vis_mon, vis_len)
vis_dat <- bind_rows(d_dep, d_mon, d_len)

# Plot!
ggplot(vis |> filter(!var == "Month"), aes(x = x, y = visregFit)) + 
  facet_grid(var ~ species) + 
  geom_point(data = vis_dat |> filter(!var == "Month"), aes(x = x, y = visregRes),
             alpha = 0.3, color = "gray50", size = 0.2) +
  geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.3, color = NA) +
  geom_line(color = "steelblue", linewidth = 1) + 
  labs(x = "Scaled variable", y = "Prediction") +
  NULL

ggsave(paste0(home, "/figures/supp/conditional.pdf"), width = 17, height = 11, units = "cm")

# Now do month (categorial)
ggplot(vis |> filter(var == "Month"), aes(x = as.factor(x), y = visregFit)) + 
  facet_wrap(~ species) + # free grid?
  geom_jitter(data = vis_dat |> filter(var == "Month"), aes(x = as.factor(x), y = visregRes),
              alpha = 0.15, color = "gray60", size = 0.2, height = 0) +
  geom_errorbar(aes(ymin = visregLwr, ymax = visregUpr), alpha = 0.6, color = "steelblue",
                width = 0, linewidth = 1) +
  geom_point(color = "steelblue", size = 2) + 
  labs(x = "Month", y = "Prediction") +
  theme(aspect.ratio = 5/6) +
  NULL

ggsave(paste0(home, "/figures/supp/conditional_month.pdf"), width = 17, height = 6, units = "cm")

Calculate indices

# Need to refit the models with the extra time argument. The reason we don't do that before is because it creates a mismatch in residual dimensions and data
# This is for interpolating between year using the random walk
extra_time <- pred_grid |> filter(!year %in% df$year) |> distinct(year) |> pull(year)

filter: removed 322,828 rows (76%), 102,718 rows remaining

distinct: removed 102,711 rows (>99%), 7 rows remaining

# Spatial random field
fit_spr_m2_extra <- sdmTMB(
  FR_spr ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  extra_time = extra_time,
  spatial = "on",
  spatiotemporal = "off",
  family = tweedie(link = "log"))

fit_her_m2_extra <- sdmTMB(
  FR_her ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  extra_time = extra_time,
  spatial = "on",
  spatiotemporal = "off",
  family = tweedie(link = "log"))

fit_sad_m2_extra <- sdmTMB(
  FR_sad ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  extra_time = extra_time,
  spatial = "on",
  spatiotemporal = "off",
  family = tweedie(link = "log"))

# Predict on grid, for indices and maps
pred_spr <- predict(fit_spr_m2_extra, newdata = pred_grid, return_tmb_object = TRUE)
pred_her <- predict(fit_her_m2_extra, newdata = pred_grid, return_tmb_object = TRUE)
pred_sad <- predict(fit_sad_m2_extra, newdata = pred_grid, return_tmb_object = TRUE)

# Make temporal index!
ncells <- filter(pred_grid, year == max(pred_grid$year)) |> nrow()

filter: removed 410,872 rows (97%), 14,674 rows remaining

index_spr <- get_index(pred_spr, area = rep(1/ncells, nrow(pred_spr$data)), bias_correct = TRUE)
index_her <- get_index(pred_her, area = rep(1/ncells, nrow(pred_her$data)), bias_correct = TRUE)
index_sad <- get_index(pred_sad, area = rep(1/ncells, nrow(pred_sad$data)), bias_correct = TRUE)

# Make long
index <- bind_rows(index_spr |> mutate(prey = "Sprat"),
                   index_her |> mutate(prey = "Herring"),
                   index_sad |> mutate(prey = "Saduria")) |> 
  mutate(Observed = ifelse(year %in% extra_time, "No", "Yes"))

mutate: new variable 'prey' (character) with one unique value and 0% NA

mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA

mutate: new variable 'Observed' (character) with 2 unique values and 0% NA

# For comparison with data
df_sum <- df |>
  pivot_longer(c("FR_spr", "FR_her", "FR_sad"), names_to = "prey", values_to = "fr") |> # because the model cannot predict the large tails!
  group_by(prey) |> 
  mutate(upr = quantile(fr, probs = 0.99)) |> 
  ungroup() |> 
  filter(fr < upr) |> 
  group_by(year, prey) |>
  summarise(mean_fr = mean(fr, na.rm = TRUE)) |> 
  mutate(prey = ifelse(prey == "FR_spr", "Sprat", prey),
         prey = ifelse(prey == "FR_her", "Herring", prey),
         prey = ifelse(prey == "FR_sad", "Saduria", prey))

pivot_longer: reorganized (FR_sad, FR_spr, FR_her) into (prey, fr) [was 12425x25, now 37275x24]

group_by: one grouping variable (prey)

mutate (grouped): new variable 'upr' (double) with 3 unique values and 0% NA

ungroup: no grouping variables

filter: removed 375 rows (1%), 36,900 rows remaining

group_by: 2 grouping variables (year, prey)

summarise: now 66 rows and 3 columns, one group variable remaining (year)

mutate (grouped): changed 66 values (100%) of 'prey' (0 new NA)

index |> 
  ggplot(aes(year, est, fill = Observed)) +
  geom_point(shape = 21, alpha = 0.7) +
  scale_fill_manual(values = c("white", "grey10")) +
  facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3, scales = "free") + 
  geom_errorbar(aes(ymin = lwr, ymax = upr), alpha = 0.4, width = 0) +
  geom_point(data = df_sum, aes(year, mean_fr, color = "Data"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
  labs(x = "Year", y = "Per capita predation", color = "") +
  scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
  theme_sleek(base_size = 9) +
  theme(legend.position = "bottom")

ggsave(paste0(home, "/figures/supp/index_ci.pdf"), width = 17, height = 7, units = "cm")

index |> 
  mutate(upr = ifelse(prey == "Sprat" & upr > 0.015, 0.015, upr)) |> 
  ggplot(aes(year, est)) +
  geom_point(alpha = 0.7) + 
  stat_smooth(method = "gam", formula = y~s(x, k=3), color = "steelblue") +
  facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3, scales = "free") + 
  labs(x = "Year", y = "Per capita predation", color = "") +
  scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
  theme_sleek(base_size = 9) +
  theme(legend.position = c(0.9, 0.9))

mutate: changed 7 values (8%) of 'upr' (0 new NA)

ggsave(paste0(home, "/figures/index.pdf"), width = 17, height = 7, units = "cm")

Sensitivity with respect to prediction grid

nd_year <- data.frame(year = unique(pred_grid$year),
                      depth_sc = 0, 
                      pred_length_sc = 0, 
                      month_f = as.factor(11))

year_pred_spr <- predict(fit_spr_m2_extra, newdata = nd_year, re_form = NA, re_form_iid = NA, se_fit = TRUE)
year_pred_her <- predict(fit_her_m2_extra, newdata = nd_year, re_form = NA, re_form_iid = NA, se_fit = TRUE)
year_pred_sad <- predict(fit_sad_m2_extra, newdata = nd_year, re_form = NA, re_form_iid = NA, se_fit = TRUE)

year_cond <- bind_rows(year_pred_spr |> mutate(prey = "Sprat"),
                       year_pred_her |> mutate(prey = "Herring"),
                       year_pred_sad |> mutate(prey = "Saduria"))

mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA

year_cond <- year_cond |> mutate(data = ifelse(year %in% extra_time, "Missing", "Present"))

mutate: new variable 'data' (character) with 2 unique values and 0% NA

# Save conditional effect of year
p1 <- ggplot(year_cond |> filter(prey == "Sprat"), aes(year, exp(est), shape = data)) +
  geom_errorbar(aes(year, ymax = exp(est + 1.96*est_se), ymin = exp(est - 1.96*est_se)), alpha = 0.2, width = 0,
                color = "steelblue3") + 
  geom_point(color = "steelblue3") + 
  #facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3, scales = "free") +
  facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3) + 
  labs(x = "", y = "", color = "") +
  scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
  scale_shape_manual(values = c(21, 16)) +
  guides(shape = "none") +
  coord_cartesian(ylim = c(0, 0.017)) +
  theme_sleek(base_size = 9) +
  theme(legend.position = c(0.9, 0.9))

filter: removed 58 rows (67%), 29 rows remaining

p1

p2 <- ggplot(year_cond |> filter(prey == "Herring"), aes(year, exp(est), shape = data)) +
  geom_errorbar(aes(year, ymax = exp(est + 1.96*est_se), ymin = exp(est - 1.96*est_se)), alpha = 0.2, width = 0,
                color = "steelblue3") + 
  geom_point(color = "steelblue3") + 
  #facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3, scales = "free") +
  facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3) + 
  labs(x = "", y = "Per capita predation", shape = "") +
  scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
  scale_shape_manual(values = c(21, 16)) +
  coord_cartesian(ylim = c(0, 0.008)) +
  theme_sleek(base_size = 9) +
  theme(legend.position = c(0.9, 0.9))

filter: removed 58 rows (67%), 29 rows remaining

p2

p3 <- ggplot(year_cond |> filter(prey == "Saduria"), aes(year, exp(est), shape = data)) +
  geom_errorbar(aes(year, ymax = exp(est + 1.96*est_se), ymin = exp(est - 1.96*est_se)), alpha = 0.2, width = 0,
                color = "steelblue3") + 
  geom_point(color = "steelblue3") + 
  #facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3, scales = "free") +
  facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3) + 
  labs(x = "Year", y = "", color = "") +
  scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
  scale_shape_manual(values = c(21, 16)) +
  guides(shape = "none") +
  coord_cartesian(ylim = c(0, 0.0012)) +
  theme_sleek(base_size = 9) +
  theme(legend.position = c(0.9, 0.9))

filter: removed 58 rows (67%), 29 rows remaining

p3

# p1 + p2 + p3
# 
# ggsave(paste0(home, "/figures/conditional_year.pdf"), width = 17, height = 7, units = "cm")

p1 / p2 / p3

ggsave(paste0(home, "/figures/conditional_year2.pdf"), width = 11, height = 17, units = "cm")




# Predict on full grid (not depth-trimmed) and calculate index
pred_grid2 <- 
  bind_rows(read_csv(paste0(home, "/data/clean/pred_grid_(1_2).csv")),
            read_csv(paste0(home, "/data/clean/pred_grid_(2_2).csv"))) |> 
  filter(quarter == 4) |> # Not needed in theory for saduria...
  mutate(depth_sc = (depth - mean(df$depth)) / sd(df$depth),
         year_f = as.factor(year),
         month_f = as.factor(3),
         year = as.integer(year),
         pred_length_sc = 0,
         ices_rect = as.factor(ices_rect)) |> 
  droplevels()

Rows: 457436 Columns: 15

── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (3): substrate, month_year, ices_rect
dbl (12): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_div

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 490110 Columns: 15
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (3): substrate, month_year, ices_rect
dbl (12): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_div

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 473,773 rows (50%), 473,773 rows remaining

mutate: converted 'year' from double to integer (0 new NA)

        converted 'ices_rect' from character to factor (0 new NA)

        new variable 'depth_sc' (double) with 7,602 unique values and 0% NA

        new variable 'year_f' (factor) with 29 unique values and 0% NA

        new variable 'month_f' (factor) with one unique value and 0% NA

        new variable 'pred_length_sc' (double) with one unique value and 0% NA

# Predict on grid, for indices and maps
pred_spr2 <- predict(fit_spr_m2_extra, newdata = pred_grid2, return_tmb_object = TRUE)
pred_her2 <- predict(fit_her_m2_extra, newdata = pred_grid2, return_tmb_object = TRUE)
pred_sad2 <- predict(fit_sad_m2_extra, newdata = pred_grid2, return_tmb_object = TRUE)

# Make temporal index!
ncells2 <- filter(pred_grid2, year == max(pred_grid2$year)) |> nrow()

filter: removed 457,436 rows (97%), 16,337 rows remaining

index_spr2 <- get_index(pred_spr2, area = rep(1/ncells2, nrow(pred_spr2$data)), bias_correct = TRUE)
index_her2 <- get_index(pred_her2, area = rep(1/ncells2, nrow(pred_her2$data)), bias_correct = TRUE)
index_sad2 <- get_index(pred_sad2, area = rep(1/ncells2, nrow(pred_sad2$data)), bias_correct = TRUE)

# Make long
index2 <- bind_rows(index_spr2 |> mutate(prey = "Sprat"),
                    index_her2 |> mutate(prey = "Herring"),
                    index_sad2 |> mutate(prey = "Saduria")) |> 
  mutate(Observed = ifelse(year %in% extra_time, "No", "Yes"))

mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'Observed' (character) with 2 unique values and 0% NA

# Calculate the index for the most common ices rectangles, the ones that cumulatively make up most stomach samples!
common_rects <- df |> 
  mutate(n = n()) |> 
  group_by(ices_rect) |> 
  summarise(n_rec = n()) |>
  ungroup() |> 
  arrange(desc(n_rec)) |> 
  mutate(cumsum = cumsum(n_rec),
         sum = sum(n_rec),
         perc = cumsum / sum) |> 
  filter(perc < 0.9)

mutate: new variable 'n' (integer) with one unique value and 0% NA
group_by: one grouping variable (ices_rect)
summarise: now 40 rows and 2 columns, ungrouped
ungroup: no grouping variables
mutate: new variable 'cumsum' (integer) with 40 unique values and 0% NA
        new variable 'sum' (integer) with one unique value and 0% NA
        new variable 'perc' (double) with 40 unique values and 0% NA
filter: removed 27 rows (68%), 13 rows remaining

common_rects

# A tibble: 13 × 5
   ices_rect n_rec cumsum   sum  perc
   <fct>     <int>  <int> <int> <dbl>
 1 42H0       3017   3017 12425 0.243
 2 40H0       1573   4590 12425 0.369
 3 43H0       1398   5988 12425 0.482
 4 43H1       1272   7260 12425 0.584
 5 41G9       1211   8471 12425 0.682
 6 41H0        630   9101 12425 0.732
 7 44H1        605   9706 12425 0.781
 8 38G5        317  10023 12425 0.807
 9 40G9        315  10338 12425 0.832
10 40G6        245  10583 12425 0.852
11 41G8        220  10803 12425 0.869
12 39G7        178  10981 12425 0.884
13 39G5        146  11127 12425 0.896

pred_grid3 <- filter(pred_grid, ices_rect %in% common_rects$ices_rect)

filter: removed 308,038 rows (72%), 117,508 rows remaining

plot_map + 
  geom_raster(data = pred_grid3 |> filter(year == 2000),
              aes(X*1000, Y*1000, fill = ices_rect)) + 
  geom_point(data = dd_plot, aes(X*1000, Y*1000, size = sample_size), alpha = 0.5) +
  scale_size(range = c(.01, 2), name = "# stomachs")

filter: removed 113,456 rows (97%), 4,052 rows remaining

Warning: Removed 2 rows containing missing values (`geom_point()`).

Warning: colourbar guide needs continuous scales.

pred_spr3 <- predict(fit_spr_m2_extra, newdata = pred_grid3, return_tmb_object = TRUE)
pred_her3 <- predict(fit_her_m2_extra, newdata = pred_grid3, return_tmb_object = TRUE)
pred_sad3 <- predict(fit_sad_m2_extra, newdata = pred_grid3, return_tmb_object = TRUE)

# Make temporal index!
ncells3 <- filter(pred_grid3, year == max(pred_grid3$year)) |> nrow()

filter: removed 113,456 rows (97%), 4,052 rows remaining

index_spr3 <- get_index(pred_spr3, area = rep(1/ncells3, nrow(pred_spr3$data)), bias_correct = TRUE)
index_her3 <- get_index(pred_her3, area = rep(1/ncells3, nrow(pred_her3$data)), bias_correct = TRUE)
index_sad3 <- get_index(pred_sad3, area = rep(1/ncells3, nrow(pred_sad3$data)), bias_correct = TRUE)

# Make long
index3 <- bind_rows(index_spr3 |> mutate(prey = "Sprat"),
                    index_her3 |> mutate(prey = "Herring"),
                    index_sad3 |> mutate(prey = "Saduria")) |> 
  mutate(Observed = ifelse(year %in% extra_time, "No", "Yes"))

mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'Observed' (character) with 2 unique values and 0% NA

index |> 
  ggplot(aes(year, est, fill = Observed)) +
  geom_point(shape = 21, alpha = 0.7) +
  scale_fill_manual(values = c("white", "grey10")) +
  facet_wrap(~prey, ncol = 1, scales = "free") + 
  geom_errorbar(aes(ymin = lwr, ymax = upr), alpha = 0.4, width = 0) +
  geom_point(data = df_sum, aes(year, mean_fr, color = "Data"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
  geom_point(data = year_cond, aes(year, exp(est), color = "Conditional Year"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
  geom_ribbon(data = year_cond, aes(year, exp(est), ymin = exp(est - 1.96*est_se), ymax = exp(est + 1.96*est_se),
                                    color = "Conditional Year"), alpha = 0.3, inherit.aes = FALSE) +
  geom_point(data = index2, aes(year, est, color = "Full grid index"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
  geom_point(data = index3, aes(year, est, color = "Common grid index"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
  labs(x = "Year", y = "Feeding ratio", color = "") +
  scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
  theme_sleek(base_size = 9) +
  theme(legend.position = "bottom")

ggsave(paste0(home, "/figures/supp/index_ci_sens.pdf"), width = 15, height = 30, units = "cm")

# Is the mismatch between data and conditional simply because data happen to be in saduria-rich areas?
saduria <- terra::rast(paste0(home, "/data/saduria-tif/FWBiomassm_raster_19812019presHighweightcor_no0_newZi.tif"))
WGS84 <- "+proj=longlat +datum=WGS84"
saduria_latlon <- terra::project(saduria, WGS84)
density_saduria <- terra::extract(saduria_latlon, pred_grid %>% dplyr::select(lon, lat), method = "bilinear")

pred_grid$density_saduria <- density_saduria$FWBiomassm_raster_19812019presHighweightcor_no0_newZi

plot_map_fc + 
  geom_raster(data = pred_grid, aes(X*1000, Y*1000, fill = density_saduria)) + 
  geom_point(data = dd_plot, aes(X*1000, Y*1000, size = sample_size), alpha = 0.5, shape = 21, color = "pink") +
  scale_size(range = c(.01, 2), name = "# stomachs") + 
  scale_fill_viridis() + 
  theme_sleek(base_size = 6) +
  facet_wrap(~year, ncol = 5)

Warning: Removed 16443 rows containing missing values (`geom_raster()`).

Warning: Removed 2 rows containing missing values (`geom_point()`).

ggsave(paste0(home, "/figures/supp/index_ci_saduria_map.pdf"), width = 17, height = 17, units = "cm")

Warning: Removed 16443 rows containing missing values (`geom_raster()`).
Removed 2 rows containing missing values (`geom_point()`).

# What is the average saduria density at the location of the hauls over year?
df_sad <- df

density_saduria <- terra::extract(saduria_latlon, df_sad %>% dplyr::select(lon, lat), method = "bilinear")

df_sad$density_saduria <- density_saduria$FWBiomassm_raster_19812019presHighweightcor_no0_newZi

df_sad <- df_sad |> 
  group_by(year) |> 
  summarise(mean_saduria = mean(density_saduria))

group_by: one grouping variable (year)
summarise: now 22 rows and 2 columns, ungrouped

p1 <- ggplot(df_sad, aes(year, mean_saduria)) + 
  geom_point() + 
  geom_line() + 
  NULL

p2 <- ggplot(year_cond |> filter(prey == "Saduria"), aes(year, exp(est))) + 
  geom_line() +
  geom_point()

filter: removed 58 rows (67%), 29 rows remaining

p1 / p2

year_cond2 <- year_cond |> filter(year %in% c(df_sad$year) & prey == "Saduria")

filter: removed 65 rows (75%), 22 rows remaining

plot(df_sad$mean_saduria ~ exp(year_cond2$est))

# Somewhat good correlation but not perfect. This is not a direct test though, because it assumes they eat saduria in proportion...

# Lastly, calculate the index for the rectangles that we have samples in each year...
# Refit without extra time, because when we filter year_rectangles that are only in data, we don't have all years and cant fit to all years...
fit_sad_m2_extra4 <- sdmTMB(
  FR_sad ~ 0 + depth_sc + s(pred_length_sc, k=3) + (1|month_f),
  data = df,
  mesh = mesh,
  time = "year",
  time_varying = ~1,
  spatial = "on",
  spatiotemporal = "off",
  family = tweedie(link = "log"))

Detected irregular time spacing with an AR(1) or random walk process.
Consider filling in the missing time slices with `extra_time`.
`extra_time = c(2010, 2011, 2015, 2016, 2017, 2019, 2020)`

summary(fit_sad_m2_extra4)

Spatial model fit by ML ['sdmTMB']
Formula: FR_sad ~ 0 + depth_sc + s(pred_length_sc, k = 3) + (1 | month_f)
Mesh: mesh (isotropic covariance)
Time column: year
Data: df
Family: tweedie(link = 'log')
 
                coef.est coef.se
depth_sc           -0.01    0.15
spred_length_sc     0.19    0.05

Smooth terms:
                    Std. Dev.
sds(pred_length_sc)         0

Random intercepts:
        Std. Dev.
month_f      0.95

Time-varying parameters:
                 coef.est coef.se
(Intercept)-1993    -9.19    1.08
(Intercept)-1994    -9.22    1.08
(Intercept)-1995    -8.63    1.01
(Intercept)-1996    -8.61    1.00
(Intercept)-1997    -8.61    1.01
(Intercept)-1998    -9.71    1.11
(Intercept)-1999    -8.58    0.99
(Intercept)-2000    -8.16    0.98
(Intercept)-2001    -7.53    0.98
(Intercept)-2002    -9.07    1.00
(Intercept)-2003    -8.43    0.98
(Intercept)-2004    -8.44    0.98
(Intercept)-2005    -8.24    0.98
(Intercept)-2006    -8.01    0.97
(Intercept)-2007    -8.80    0.98
(Intercept)-2008    -8.67    0.98
(Intercept)-2009    -8.61    0.99
(Intercept)-2012    -7.96    0.99
(Intercept)-2013    -8.48    0.99
(Intercept)-2014    -8.22    0.98
(Intercept)-2018    -7.54    0.95
(Intercept)-2021    -7.14    0.94

Dispersion parameter: 0.73
Tweedie p: 1.54
Matérn range: 80.79
Spatial SD: 2.99
ML criterion at convergence: -782.798

See ?tidy.sdmTMB to extract these values as a data frame.

**Possible issues detected! Check output of sanity().**

df_saduria <- df |> mutate(year_rect_id = paste(year, ices_rect, sep = "."))

mutate: new variable 'year_rect_id' (character) with 159 unique values and 0% NA

pred_grid4 <- pred_grid |> 
  mutate(year_rect_id = paste(year, ices_rect, sep = ".")) |> 
  filter(year_rect_id %in% c(df_saduria$year_rect_id))

mutate: new variable 'year_rect_id' (character) with 1,769 unique values and 0% NA
filter: removed 383,605 rows (90%), 41,941 rows remaining

plot_map_fc + 
  geom_raster(data = pred_grid4, aes(X*1000, Y*1000, fill = ices_rect)) + 
  facet_wrap(~year) + 
  geom_point(data = df, aes(X*1000, Y*1000))

Warning: Removed 1249 rows containing missing values (`geom_raster()`).

Warning: Removed 39 rows containing missing values (`geom_point()`).

Warning: colourbar guide needs continuous scales.

pred_sad4 <- predict(fit_sad_m2_extra4, newdata = pred_grid4, return_tmb_object = TRUE)

ncells4 <- filter(pred_grid4, year == max(pred_grid4$year)) |> nrow()

filter: removed 33,599 rows (80%), 8,342 rows remaining

index_sad4 <- get_index(pred_sad4, area = rep(1/ncells4, nrow(pred_sad4$data)), bias_correct = TRUE)

manual_index <- pred_sad4$data |> 
  group_by(year) |> 
  summarise(index_manual = mean(exp(est)))

group_by: one grouping variable (year)
summarise: now 22 rows and 2 columns, ungrouped

index_sad4 |> 
  ggplot(aes(year, est)) +
  geom_point() +
  geom_errorbar(aes(ymin = lwr, ymax = upr), alpha = 0.4, width = 0) +
  geom_point(data = df_sum |> filter(prey == "Saduria"), aes(year, mean_fr, color = "Data"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
  geom_point(data = manual_index, aes(year, index_manual, color = "Manual index"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
  geom_point(data = year_cond |> filter(prey == "Saduria"), aes(year, exp(est), color = "Conditional Year"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
  geom_ribbon(data = year_cond |> filter(prey == "Saduria"), aes(year, exp(est), ymin = exp(est - 1.96*est_se), ymax = exp(est + 1.96*est_se),
                                    color = "Conditional Year"), alpha = 0.3, inherit.aes = FALSE) +
  labs(x = "Year", y = "Feeding ratio", color = "") +
  scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
  coord_cartesian(ylim = c(0, 0.015)) +
  theme(legend.position = "bottom")

filter (grouped): removed 44 rows (67%), 22 rows remaining
filter: removed 58 rows (67%), 29 rows remaining
filter: removed 58 rows (67%), 29 rows remaining

# The only thing that could change the conditional from the index is depth I guess? And size
t <- pred_grid4 |> 
  group_by(year) |> 
  summarise(mean_depth_sc = mean(depth_sc))

group_by: one grouping variable (year)
summarise: now 22 rows and 2 columns, ungrouped

df |> 
  group_by(year) |> 
  summarise(mean = mean(pred_length_sc)) |> 
  ggplot(aes(year, mean)) + 
  geom_point()

group_by: one grouping variable (year)
summarise: now 22 rows and 2 columns, ungrouped

# Make a new conditional with values closer to the prediction grid
nd_year2 <- data.frame(year = unique(pred_grid4$year),
                       depth_sc = t$mean_depth_sc, 
                       pred_length_sc = 0, 
                       month_f = as.factor(11))

year_pred_sad2 <- predict(fit_sad_m2_extra4, newdata = nd_year2, re_form = NA, re_form_iid = NULL, se_fit = TRUE)

index_sad4 |> 
  ggplot(aes(year, est)) +
  geom_point() +
  geom_errorbar(aes(ymin = lwr, ymax = upr), alpha = 0.4, width = 0) +
  geom_point(data = df_sum |> filter(prey == "Saduria"), aes(year, mean_fr, color = "Data"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
  geom_point(data = manual_index, aes(year, index_manual, color = "Manual index"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
  geom_point(data = year_pred_sad2, aes(year, exp(est), color = "Conditional Year"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
  geom_ribbon(data = year_pred_sad2, aes(year, exp(est), ymin = exp(est - 1.96*est_se), ymax = exp(est + 1.96*est_se),
                                    color = "Conditional Year"), alpha = 0.3, inherit.aes = FALSE) +
  labs(x = "Year", y = "Feeding ratio", color = "") +
  scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
  theme(legend.position = "bottom")

filter (grouped): removed 44 rows (67%), 22 rows remaining

index_sad4 |> 
  ggplot(aes(year, est)) +
  geom_line() +
  geom_line(data = year_pred_sad2, aes(year, exp(est), color = "Conditional Year")) +
  geom_line(data = year_pred_sad, aes(year, exp(est), color = "Conditional Year fixed")) +
  labs(x = "Year", y = "Feeding ratio", color = "") +
  scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
  theme(legend.position = "bottom")

# Ok, now I can get the index and the conditional effect match, i.e., by setting depths to match. Can I get conditional to match data by setting depth, month and size to match?
tt <- df |> 
  group_by(month_f, year) |> 
  summarise(n = n()) |> 
  group_by(year) |> 
  filter(n == max(n)) |> 
  arrange(year)

group_by: 2 grouping variables (month_f, year)
summarise: now 43 rows and 3 columns, one group variable remaining (month_f)
group_by: one grouping variable (year)
filter (grouped): removed 21 rows (49%), 22 rows remaining

t <- df |> 
  group_by(year) |> 
  summarise(mean_depth_sc = mean(depth_sc),
            mean_pred_length_sc = mean(pred_length_sc))

group_by: one grouping variable (year)
summarise: now 22 rows and 3 columns, ungrouped

nd_year3 <- data.frame(year = unique(df$year),
                       depth_sc = t$mean_depth_sc, 
                       pred_length_sc = t$mean_pred_length_sc, 
                       month_f = tt$month_f) # Most common month!

year_pred_sad3 <- predict(fit_sad_m2_extra4, newdata = nd_year3, re_form = NA, re_form_iid = NULL, se_fit = TRUE)

ggplot() +
  geom_point(data = df_sum |> filter(prey == "Saduria"), aes(year, mean_fr, color = "Data"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
  geom_line(data = year_pred_sad3, aes(year, exp(est), color = "Conditional Year"), alpha = 0.6, inherit.aes = FALSE, shape = 2) +
  geom_ribbon(data = year_pred_sad3, aes(year, exp(est), ymin = exp(est - 1.96*est_se), ymax = exp(est + 1.96*est_se),
                                         color = "Conditional Year"), alpha = 0.3, inherit.aes = FALSE) +
  labs(x = "Year", y = "Feeding ratio", color = "") +
  scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
  theme(legend.position = "bottom")

filter (grouped): removed 44 rows (67%), 22 rows remaining

Warning in geom_line(data = year_pred_sad3, aes(year, exp(est), color =
"Conditional Year"), : Ignoring unknown parameters: `shape`

# Is the mean higher because I have tails that the model doesn't capture?
df_sum2 <- df |>
  filter(FR_sad < quantile(FR_sad, prob = 0.95)) |> 
  group_by(year) |>
  summarise(Saduria = mean(FR_sad))

filter: removed 622 rows (5%), 11,803 rows remaining
group_by: one grouping variable (year)
summarise: now 22 rows and 2 columns, ungrouped

ggplot() +
  geom_line(data = year_pred_sad3, aes(year, exp(est), color = "Conditional Year"), alpha = 0.6, inherit.aes = FALSE) +
  geom_ribbon(data = year_pred_sad3, aes(year, exp(est), ymin = exp(est - 1.96*est_se), ymax = exp(est + 1.96*est_se),
                                    fill = "Conditional Year"), alpha = 0.3, inherit.aes = FALSE, color = NA) +
  geom_line(data = df_sum |> filter(prey == "Saduria"), aes(year, mean_fr, color = "Data"), alpha = 0.6, inherit.aes = FALSE) +
  geom_line(data = df_sum2, aes(year, Saduria, color = "Data trimmed"), alpha = 0.6, inherit.aes = FALSE) +
  labs(x = "Year", y = "Feeding ratio", color = "") +
  scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
  theme(legend.position = "bottom")

filter (grouped): removed 44 rows (67%), 22 rows remaining

# YEEEESS FINALLY!!!
# TODO: still I want to estimate the effect of index prediction vs conditional... 

Plot spatial predictions

spatial_preds <- bind_rows(pred_spr$data |> mutate(prey = "Sprat"),
                           pred_her$data |> mutate(prey = "Herring"),
                           pred_sad$data |> mutate(prey = "Saduria"))

mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA
mutate: new variable 'prey' (character) with one unique value and 0% NA

# Plot spatial random effect
plot_map +
  geom_raster(data = spatial_preds |> filter(year == 2000),
              aes(X*1000, Y*1000, fill = omega_s)) +
  scale_fill_gradient2(name = "Spatial random field") +
  facet_wrap(~prey) +
  geom_sf(size = 0.1) + 
  theme_sleek(base_size = 9) + 
  theme(legend.position = c(0.1, 0.7)) +
  theme(legend.key.height = unit(0.6, "line"),
        legend.key.width = unit(0.2, "line"))

filter: removed 1,232,616 rows (97%), 44,022 rows remaining

Warning: Removed 1701 rows containing missing values (`geom_raster()`).

ggsave(paste0(home, "/figures/supp/diet_omega_s.pdf"), width = 17, height = 6, units = "cm")

Warning: Removed 1701 rows containing missing values (`geom_raster()`).

# Plot spatiotemporal predictions
plot_map_fc +
  geom_raster(data = spatial_preds |> filter(prey == "Sprat"), aes(X*1000, Y*1000, fill = exp(est))) +
  scale_fill_viridis(trans = "sqrt", name = "FR sprat") +
  facet_wrap(~year) +
  geom_sf(size = 0.1)

filter: removed 851,092 rows (67%), 425,546 rows remaining

Warning: Removed 16443 rows containing missing values (`geom_raster()`).

ggsave(paste0(home, "/figures/supp/sprat_spatiotemporal_diet.pdf"), width = 17, height = 17, units = "cm")

Warning: Removed 16443 rows containing missing values (`geom_raster()`).

plot_map_fc +
  geom_raster(data = spatial_preds |> filter(prey == "Herring"), aes(X*1000, Y*1000, fill = exp(est))) +
  scale_fill_viridis(trans = "sqrt", name = "FR herring") +
  facet_wrap(~year) +
  geom_sf(size = 0.1)

filter: removed 851,092 rows (67%), 425,546 rows remaining

Warning: Removed 16443 rows containing missing values (`geom_raster()`).

ggsave(paste0(home, "/figures/supp/herring_spatiotemporal_diet.pdf"), width = 17, height = 17, units = "cm")

Warning: Removed 16443 rows containing missing values (`geom_raster()`).

plot_map_fc +
  geom_raster(data = spatial_preds |> filter(prey == "Saduria"), aes(X*1000, Y*1000, fill = exp(est))) +
  scale_fill_viridis(trans = "sqrt", name = "FR saduria") +
  facet_wrap(~year) +
  geom_sf(size = 0.1)

filter: removed 851,092 rows (67%), 425,546 rows remaining

Warning: Removed 16443 rows containing missing values (`geom_raster()`).

ggsave(paste0(home, "/figures/supp/saduria_spatiotemporal_diet.pdf"), width = 17, height = 17, units = "cm")

Warning: Removed 16443 rows containing missing values (`geom_raster()`).

# Plot spatiotemporal predictions for a year and all species
plot_map +
  geom_raster(data = spatial_preds|> filter(year == "2000"), aes(X*1000, Y*1000, fill = exp(est))) +
  scale_fill_viridis(trans = "sqrt", name = "Per cap. predation") +
  facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria"))) +
  geom_sf(size = 0.1) + 
  theme_sleek(base_size = 9) + 
  theme(legend.position = c(0.08, 0.7),
        legend.key.height = unit(0.6, "line"),
        legend.key.width = unit(0.2, "line"))

filter: removed 1,232,616 rows (97%), 44,022 rows remaining

Warning: Removed 1701 rows containing missing values (`geom_raster()`).

ggsave(paste0(home, "/figures/spatial_prey_prediction_2000.pdf"), width = 17, height = 6, units = "cm")

Warning: Removed 1701 rows containing missing values (`geom_raster()`).

Plot spatial overlap vs and populationa and per capita feeding

# Join feeding ratio indices (per capita) and overlap
cod_spr <- read_csv(paste0(home, "/output/cod_pel_sum_ovrlap.csv")) |>
  dplyr::select(year, cod_spr_ovr_tot) |> 
  rename(overlap = cod_spr_ovr_tot) |> 
  mutate(prey = "Sprat")

Rows: 27 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (5): year, cod_spr_ovr_tot, cod_spr_ovr_sd_tot, cod_her_ovr_tot, cod_her...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (overlap)

mutate: new variable 'prey' (character) with one unique value and 0% NA

cod_her <- read_csv(paste0(home, "/output/cod_pel_sum_ovrlap.csv")) |>
  dplyr::select(year, cod_her_ovr_tot) |> 
  rename(overlap = cod_her_ovr_tot) |> 
  mutate(prey = "Herring")

Rows: 27 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (5): year, cod_spr_ovr_tot, cod_spr_ovr_sd_tot, cod_her_ovr_tot, cod_her...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (overlap)

mutate: new variable 'prey' (character) with one unique value and 0% NA

cod_ben <- read_csv(paste0(home, "/output/cod_ben_sum_ovrlap.csv")) |>
  dplyr::select(year, cod_sad_ovr_tot) |> 
  rename(overlap = cod_sad_ovr_tot) |> 
  mutate(prey = "Saduria")

Rows: 54 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (3): year, quarter, cod_sad_ovr_tot

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (overlap)

mutate: new variable 'prey' (character) with one unique value and 0% NA

overlap <- bind_rows(cod_spr, cod_her, cod_ben)

index_ovr <- index |>
  left_join(overlap, by = c("year", "prey")) |> 
  drop_na()

left_join: added one column (overlap)
           > rows only in x     6
           > rows only in y  (  0)
           > matched rows     108    (includes duplicates)
           >                 =====
           > rows total       114
drop_na: removed 6 rows (5%), 108 rows remaining

# Join in cod biomass data to calculate snapshot predation
cod_pred <- read_csv(paste0(home, "/output/pred_cod.csv")) |> 
  filter(quarter == 4)

Rows: 943428 Columns: 38
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (3): substrate, month_year, ices_rect
dbl (35): X, Y, year, lon, lat, depth, quarter, month, oxy, temp, sal, sub_d...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 471,714 rows (50%), 471,714 rows remaining

# Left join cod biomass density onto the grid of diet predictions
spatial_preds_dens <- left_join(spatial_preds,
                                cod_pred |> dplyr::select(est, X, Y, year) |> rename(est_codbiom = est))

rename: renamed one variable (est_codbiom)
Joining with `by = join_by(X, Y, year)`left_join: added one column (est_codbiom)
           > rows only in x       6,177
           > rows only in y  (   48,227)
           > matched rows     1,270,461
           >                 ===========
           > rows total       1,276,638

# Multiply local feeding ratio with biomass density of cod
spatial_preds_dens <- spatial_preds_dens |> 
  drop_na(est_codbiom) |> 
  mutate(pred = exp(est) * est_codbiom)

drop_na: removed 6,177 rows (<1%), 1,270,461 rows remaining
mutate: new variable 'pred' (double) with 1,270,461 unique values and 0% NA

# Plot map!
plot_map +
  geom_raster(data = spatial_preds_dens |> filter(year == "2000"), aes(X*1000, Y*1000, fill = pred)) +
  scale_fill_viridis(trans = "sqrt", name = "Pop. predation") +
  facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria"))) +
  geom_sf(size = 0.1) + 
  theme_sleek(base_size = 9) + 
  theme(legend.position = c(0.08, 0.7),
        legend.key.height = unit(0.6, "line"),
        legend.key.width = unit(0.2, "line"))

filter: removed 1,226,652 rows (97%), 43,809 rows remaining

Warning: Removed 1698 rows containing missing values (`geom_raster()`).

ggsave(paste0(home, "/figures/spatial_pop_predation_prediction_2000.pdf"), width = 17, height = 6, units = "cm")

Warning: Removed 1698 rows containing missing values (`geom_raster()`).

# index? Just take the average
spatial_preds_dens |> 
  group_by(year, prey) |> 
  summarise(mean_pop_pred = mean(pred)) |> 
  ggplot(aes(year, mean_pop_pred)) +
  geom_point(alpha = 0.7) + 
  stat_smooth(method = "gam", formula = y~s(x, k=3), color = "steelblue") +
  facet_wrap(~factor(prey, levels = c("Sprat", "Herring", "Saduria")), ncol = 3, scales = "free") + 
  labs(x = "Year", y = "Population predation", color = "") +
  scale_x_continuous(breaks = seq_range(unique(df$year), by = 8)) +
  theme_sleek(base_size = 9) +
  theme(legend.position = c(0.9, 0.9))

group_by: 2 grouping variables (year, prey)
summarise: now 87 rows and 3 columns, one group variable remaining (year)

ggsave(paste0(home, "/figures/pop_pred_index.pdf"), width = 17, height = 7, units = "cm")


# Summarize across years
spatial_preds_dens_sum <- spatial_preds_dens |> 
  group_by(prey, year) |> 
  summarise(tot_pred = sum(pred))

group_by: 2 grouping variables (prey, year)
summarise: now 87 rows and 3 columns, one group variable remaining (prey)

# Left_join with index data
index_ovr <- index_ovr |>
  left_join(spatial_preds_dens_sum, by = c("year", "prey")) |> 
  drop_na()

left_join: added one column (tot_pred)
           > rows only in x     0
           > rows only in y  (  6)
           > matched rows     108
           >                 =====
           > rows total       108
drop_na: no rows removed

# Check some regressions...
# summary(lm(tot_pred ~ overlap, data = index_ovr |> filter(prey == "Sprat")))
# summary(lm(tot_pred ~ overlap, data = index_ovr |> filter(prey == "Herring")))
# summary(lm(tot_pred ~ overlap, data = index_ovr |> filter(prey == "Saduria")))
# 
# summary(lm(est ~ overlap, data = index_ovr |> filter(prey == "Sprat")))
# summary(lm(est ~ overlap, data = index_ovr |> filter(prey == "Herring")))
# summary(lm(est ~ overlap, data = index_ovr |> filter(prey == "Saduria")))

# Calculate correlations between FR, predation and overlap
cor <- plyr::ddply(index_ovr, c("prey"),
                   summarise,
                   cor_fr_ovr = round(cor(overlap, est), 2),
                   cor_pred_ovr = round(cor(overlap, tot_pred), 2)) |> 
  pivot_longer(c("cor_fr_ovr", "cor_pred_ovr")) |> 
  mutate(name = ifelse(name == "cor_fr_ovr", "Per capita", "Population"))

summarise: now one row and 2 columns, ungrouped
summarise: now one row and 2 columns, ungrouped
summarise: now one row and 2 columns, ungrouped
pivot_longer: reorganized (cor_fr_ovr, cor_pred_ovr) into (name, value) [was 3x3, now 6x3]
mutate: changed 6 values (100%) of 'name' (0 new NA)

p <- index_ovr |>
  rename(fr = est) |> 
  pivot_longer(c("fr", "tot_pred")) |> 
  mutate(id = paste(name, prey, sep = ";")) %>%
  split(.$id) |>
  purrr::map(~lm(value ~ overlap, data = .x)) |>
  purrr::map_df(broom::tidy, .id = 'id') |>
  filter(term == 'overlap') |>
  separate(id, sep = ";", into = c("name", "prey"), remove = FALSE) |> 
  mutate(name = ifelse(name == "fr", "Per capita", "Population"))

rename: renamed one variable (fr)
pivot_longer: reorganized (fr, tot_pred) into (name, value) [was 108x10, now 216x10]
mutate: new variable 'id' (character) with 6 unique values and 0% NA
filter: removed 6 rows (50%), 6 rows remaining
mutate: changed 6 values (100%) of 'name' (0 new NA)

p

# A tibble: 6 × 8
  id               name       prey    term  estimate std.error statistic p.value
  <chr>            <chr>      <chr>   <chr>    <dbl>     <dbl>     <dbl>   <dbl>
1 fr;Herring       Per capita Herring over…  3.06e-4   3.11e-3   0.0982   0.923 
2 fr;Saduria       Per capita Saduria over… -2.48e-4   4.88e-4  -0.509    0.613 
3 fr;Sprat         Per capita Sprat   over…  9.58e-3   4.88e-3   1.96     0.0610
4 tot_pred;Herring Population Herring over… -2.83e+1   1.39e+4  -0.00203  0.998 
5 tot_pred;Saduria Population Saduria over…  1.31e+3   5.64e+2   2.32     0.0243
6 tot_pred;Sprat   Population Sprat   over…  2.62e+4   1.94e+4   1.35     0.188 

# Plot!
index_ovr |> 
  rename('Per capita' = est, 
         'Population' = tot_pred) |> 
  pivot_longer(c('Per capita', 'Population')) |> 
  ggplot(aes(overlap, value)) +
  geom_point(alpha = 0.5) +
  ggh4x::facet_grid2(name ~ prey, scales = "free", independent = "y") +
  # geom_text(data = cor,
  #           aes(label = paste("r=", value, sep = "")), x = -Inf, y = Inf, hjust = -.1, vjust = 2.5,
  #           inherit.aes = FALSE, fontface = "italic", size = 2.5, color = "tomato3") +
  geom_text(data = p,
            aes(label = paste("p=", round(p.value, digits = 3), sep = "")), x = -Inf, y = Inf, hjust = -.1, vjust = 2.5,
            inherit.aes = FALSE, fontface = "italic", size = 2.5, color = "tomato3") +
  labs(y = "Predation", x = "Spatial overlap", color = "Year") +
  theme_sleek(base_size = 9) +
  theme(legend.position = "bottom", aspect.ratio = 1) +
  NULL

rename: renamed 2 variables (Per capita, Population)
pivot_longer: reorganized (Per capita, Population) into (name, value) [was 108x10, now 216x10]

# Summarize the differences here, make fake facet grid free axis

ggsave(paste0(home, "/figures/fr_overlap_cor.pdf"), width = 17, height = 11, units = "cm")